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COMPUTATION  OF  LOW  SPEED  VISCOUS  FLOWS  WITH  HEAT  ADDITION 

Ashvln  Hosangadl*  and  Charles  L.  Merkle 
The  Pennsylvania  State  University 
Department  of  Mechanical  Engineering 
University  Park,  PA  16602 


The  use  of  Implicit  time  dependent  schemes  for 
the  numerical  solution  of  low  speed,  low  Reynolds 
number  flows  with  heat  addition  Is  investigated. 
Stability  analyses  show  that  the  errors  Introduced  by 
approximate  factorization  give  rise  to  Instability  at 
Reynolds  numbers  around  100.  Specifically,  It  la  the 
cross-derivative  errors  between  the  viscous  and 
lnvlscld  terms  that  cause  problems.  When  exact 
Inversion  techniques  are  used,  the  system  becomes 
strongly  stable  and  numerical  experiments  show  rapid 
convergence.  Comparisons  of  outflow  boundary 
conditions  show  that  viscous  and  lnviscld 
formulations  give  identical  results  over  a  wide  rang* 
of  Reynolds  numbers  when  buoyancy  Is  omitted,  but 
with  buoyancy  present  the  lnviscld  boundary 
conditions  are  unstable.  Flowfleld  results  for  a 
range  of  low  Reynolds  conditions  with  and  without 
buoyancy  are  given  to  show  the  manner  In  which  the 
flowfleld  changes  as  these  physical  parameters  are 
varied. 

INTRODUCTION 

Time  dependent  procedures  have  come  to  be  the 
accepted  method  for  the  numerical  solution  of 
compressible  flows.  The  strengths  of  these  methods 
arise  because  they  are  easily  expressed  In 
generalized  coordinates,  they  allow  a  single 
procedure  to  be  used  Tor  either  viscous  or  lnviscld 
flows,  and  they  allow  the  use  of  central  differences 
for  the  convective  terms  even. when  the  lnviscld 
equations  are  considered.  Their  capabilities  in 
lnvlscld  flows  make  them  particularly  suitable  for 
the  Important  class  or  high  Reynolds  number  viscous 
Hows.  Recent  publications*-4  have  reported 
effective  methods  for  extending  these  compressible 
flow  techniques  to  low  Mach  number  conditions  where 
eigenvalue  stiffness  had  heretofore  precluded 
convergence.  The  present  paper  looks  at  the 
additional  extension  to  low  Reynolds  number  flows. 

The  physical  problem  used  to  demonstrate  the  low 
Reynolds  number  conditions  Is  that  of  strong  heat 
addition  In  a  constant  area  duct.  The  source  of  this 
problem  arises  because  of  our  Interest  In  the 
absorption  of  a  laser  beam  by  a  riowlng  gas5. 

The  Impaired  convergence  of  time  dependent 
schemes  at  low  Mach  number  Is  known  to  arise  because 
of  the  stirfness  or  the  eigenvalues  and  the  singular 
behavior  of  the  pressure  gradient.  Preconditioning 
of  the  time  derivatives1-?  has  been  shown  to  make  the 
convergence  rate  Independent  of  Mach  number,  but  a 
rescaling  of  the  pressure  Is  necessary4  to  ensure 
convergence  to  tight  tolerances  at  Mach  numbers  below 
about  10-4.  This  rescaling  Is  best  accomplished  by  a 
perturbation  expansion  of  the  equations  of  motion4’5. 
The  presence  of  the  source  term  due  to  buoyancy  has 
also  been  shown^  to  require  special  treatment  when 
Froude  numbers  are  small. 


Graduate  Assistant;  Professor 


AFOSR.TR ■>  8  8  -  0  4  6  8 

The  characteristics  of  time  dependent  scnemes  at 
low  Reynolds  numbers  have  received  less  attention. 
Although  the  diffusion  provided  by  the  viscous  terms 
should  smooth  the  solution  and  provide  enhanced 
convergence,  their  presence  can  slow  convergence  or 
even  lead  to  divergence  when  approximate 
decompositions  are  used  In  conjunction  with  implicit 
formulations.  (The  characteristics  of  approximately 
factored®’®  methods  are  discussed  herein  a  l  compared 
with  results  from  direct  Inversion  procedures.)  In 
addition,  the  effects  of  Inflow  and  outflow  boundary 
conditions  are  studied  to  determine  their  effect  on 
convergence  and  solution  accuracy  as  the  Reynolds 
number  is  decreased. 

PROBLEM  FORMULATION  AND  EQUATIONS  OF  MOTION 

As  a  test  problem,  we  uoe  tne  flow  through  a 
vertical  constant  area  duct  with  specified  volumetric 
heat  addition  to  simulate  laser  absorption  in  a 
flowing  gas5.  To  parameterize  the  problem,  the  heat 
source,  Qh,  is  chosen  as  the  algebraic  function, 

Oh  -  P  Min  |2f 3/2  -  3f*i,  i|  (1) 


f(x.y)  -  |(x-xa)2  *  (y-ya)2)/|(x-xb)2  ♦  (y-yb)2]  (2) 

Here  P  is  a  scaling  parameter  indicative  of  the 
maximum  heating  rate  and  (xa,  ya)  and  (xb,  yb) 
speciry  the  coordinates  about  which  the  symmetric 
heat  addition  is  centered  and  Its  rate  or  decrease. 
For  the  present  calculations,  only  the  peak  heating 
rate  P  was  varied.  Cross-sectional  contours  or  the 
constant  heat  addition  curves  along  with  the  geometry 
or  the  problem  Is  given  In  Fig.  1. 

The  equations  used  to  describe  the  very  low 
speed  riows  of  Interest  here  are  obtained  by  the  same 
Mach  number  expansion  used  In  Ref.  A.  Extension  of 
this  expansion  to  the  viscous  terms  shows  that  the 
viscous  dissipation  function  In  the  energy  equation 
may  be  neglected  leaving  only  heat  conduction  terms 
on  the  right-hand  side  and  leads  to  no 
simplifications  In  the  right-hand  side  of  the 
momentum  equations.  The  reduced  equations  In 
generalized  coordinates  become: 

3.Q  *  3fE  *  -  H  ‘  I  VRrr3r  *  3  > 

t  C  n  C  n 

4  3n(Rnf3f  4  Rnn3jlJQ  (3) 

n  c  nn  o 

The  notation  used  for  the  lnviscld  terms  Is  the  same 
as  that  of  Ref.  A  while  the  viscous  notation  Is 
similar  to  that  of  Ref.  10  with  the  above-mentioned 
simplification.  The  source  term  H  contains  both  the 
heat  source  term  and  the  buoyancy  term. 

To  solve  Eon.  3  we  use  Euler  Implicit  time 
differencing  for  all  terms  except  for  the 
cross-derivatives  on  the  right-hand  sloe,  which  are 
differenced  explicitly.  Central  differencing  In 
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space  Is  used  Tor  both  the  convective  and  diffusive 
terns.  The  Implicit  formulation  gives  rise  to  a 
large  sparse,  pentadiagonal  matrix  whose  solution  was 
attempted  by  the  approximate  factorization  technique 
as  well  as  by  an  exact,  direct  Inversion  procedure. 

Approximate  Factorization 

The  most  common  technique  for  solving  Implicit 
algorithms  Is  the  approximate  factorization  procedure 
of  Douglas  and  Cunn^i.  In  this  approach,  the 
pentadiagonal  matrix  Is  approximated  by  two 
trldlagonal  matrices  as. 


(S  ♦ 


ata^A 


AtS^R^JlS 


-1 


IS  *  4t3  8  -  4t3  R  3  J)AQ 
n  n  nn  n 


-At  R 


where  R  Is  the  residual  of  the  operator  In  Eqn.  3 
evaluated  at  the  old  time  level,  and  S  -  I-Dit  where 
D  is  the  Jacobian  of  H. 

The  approximate  factorization  la  equivalent  to 
solving  the  original  differential  operator  In  Eqn.  3 
with  an  additive  error  wAF  given  by. 
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Note  this  error  Includes  terms  of  three  kinds.  There 
are  lnvlscld-lnvlscld  operator  products, 
viscous-viscous  operator  products  and 
viseous-lnvlscid  operator  products.  The 
lnvlscld-lnvlscld  terms  are  the  familiar  ones  from 
the  Euler  equations.  These  are  well  known  to  slow 
convergence  at  hlgn  CFL's,  but  never  cause 
Instability  or  divergence.  Similar  comments  hold  for 
the  viscous-viscous  error  term.  In  a  pure  diffusion 
equation  (no  convection)  these  terms  will  alow 
convergence  but  the  approximately  factored  scheme 
remains  unconditionally  stable. 

The  existence  of  an  unconditionally  stable 
factorization  for  the  lnvlscld  terms  alone  and  for 
the  diffusion  terms  alone  does  not,  however, 
guarantee  unconditional  stability  for  the  combined 
vlseous-lnvlscld  system  because  of  the 
vlseous-lnvlscld  error  terms  In  Eqn.  5.  Stability 
analyses  of  the  complete  vector  equations  shows  that 
the  vlseous-lnvlscld  error  terms  can  be  strongly 
destabilizing  at  moderate  (order  one)  Reynolds 
numbers.  At  high  or  very  low  Reynolds  numbers,  these 
terms  become  negligible  and  have  no  effect.  For  our 
specific  problem,  both  stability  analyses  and 
numerical  experiments  showed  that  approximate 
factorization  was  unstable  at  Reynolds  numbers  below 
500  when  the  stirfness  was  removed  from  the  lnvlscld 
terms.  Making  the  lnvlscld  operators  stiff  deferred 
this  Instability  to  lower  Reynolds  number,  but  the 
stiffness  then  caused  unacceptably  slow  convergence. 
Thus,  approximate  factorization  of  the  low  Mach 
number  equations  falls  at  low  Reynolds  numbers. 

Direct  Solution  Procedures 


Analysis  of  the  AF  scheme  reveals  Its  limited 
usefulness  at  low  Reynolds  numbers.  To  understand 
the  limitations  or  the  approximate  procedure  and  to 
provide  solutions  to  the  problem  of  Interest,  we  also 
considered  direct  solutions  of  the  operator  In  Eqn.  3- 
For  these  solutions,  the  sparse  pentadiagonal  matrix 
was  partitioned  as  a  block  trldlagonal  system  which 
was  then  solved  by  the  matrix  version  of  the  standard 


Thomas  algorithm.  The  result  requires  considerable 
machine  storage  but  provides  an  extremely  robust 
solver  as  Is  noted  below.  For  a  two-dimensional 
problem,  both  storage  and  CPU  time  are  primarily 
determined  by  the  number  of  grid  points  in  the 
smaller  direction.  Experimental  cheeks  on  a  scalar 
processing  machine  snowed  that  when  20  to  30  grid 
points  were  used  In  one  direction,  the  direct 
procedure  gave  converged  solutions  In  about  the  same 
time  as  approximate  factorization,  albeit,  at  vastly 
different  convergence  rates.  Thus,  a  erld  of  50  x  25 
will  require  nominally  the  same  time  for  the  direct 
as  for  the  AF  procedure.  Mors  advanced  sparse  macrlx 
solution  procedures12  are  available  which  would 
reduce  storage  as  compared  with  the  present 
trldlagonal  procedure.  These  should  also  reduce  the 
CPU  time  requirement. 

The  advantage  of  tne  direct  procedure  Is  Its 
Improved  robustness  and  Its  rapid  rate  of  convergence. 
Its  robustness  Is  demonstrated  by  the  fact  that 
converged  solutions  are  obtained  for  the  present 
problem,  while  appropriate  procedures  fall.  Its 
rapid  rate  of  convergence  Is  most  readily 
demonstrated  by  example.  Figure  2  snows  the 
convergence  rate  of  tne  direct  procedure  for  the 
present  problem  for  a  series  of  different  CFL's.  At 
a  CFL  of  5.  the  direct  proeecure  converges  at  about 
the  same  rate  as  an  approximate  factorization  scheme. 
(The  optimum  CFL  for  AF  schemes  Is  around  5.)  As  CFL 
Is  Increased,  the  convergence  rate  accelerates 
rapidly.  At  a  CFL  of  50,  convergence  to  machine 
accuracy  requires  only  17  Iterations  and  at  5000,  It 
Is  further  reduced  to  11  time  steps.  Tests  with 
other  grid  sizes  show  that  the  rate  of  convergence  Is 
Independent  of  the  number  or  grids  or  even  of  grid 
stretching  although  the  time  per  Iteration  Is 
strongly  related  to  the  number  or  grids.  Finally, 
comparison  of  the  present  rates  or  convergence  for  a 
two-dimensional  problem  with  those  Tor  a 
one-dlmenslonal  problem  show  that  the  convergence 
rate  is  independent  of  the  dimensionality  of  the 
problem  when  direct  Inversion  Is  used.  The  normally 
encountered  slow-down  In  convergence  rate  In 
multidimensions  as  compared  to  one  dimension  arises 
because  or  approximate  factorization  errors. 

Boundary  Conditions 

The  required  boundary  conditions  for  the  present 
problem  Include  the  specification  of  conditions  on 
the  wall,  on  the  axis  of  symmetry  and  on  the  upstream 
and  downstream  boundaries.  The  boundary  conditions 
on  the  wall  are  the  no-sllp  and  adiabatic  conditions. 
These  three  conditions  are  supplemented  by  an 
application  of  the  normal  momentum  equation  on  the 
wall  which  Is  applied  by  using  one-sided  differences 
for  the  normal  derivatives.  At  the  center  line, 
symmetry  conditions  are  applied.  Because  these 
conditions  are  simple  and  their  application  Is 
straightforward,  the  details  are  omitted. 

The  boundary  conditions  at  the  Inlet  and  exit 
are  of  more  Interest.  At  high  Reynolds  numbers  where 
the  flow  Is  predominantly  lnvlscld,  the  use  of 
lnvlscld  boundary  procedures  appears  appropriate.  As 
the  Reynolds  number  Is  decreased  and  viscous  effects 
begin  to  dominate,  the  arguments  for  viscous  boundary 
conditions  become  more  compelling.  The  range  of 
demarcation  between  these  two  types  of  boundary 
conditions  Is  not  Immediately  apparent  aid 
conclusions  are  based  on  hindsight  gained  from 
numerical  experimentation. 

For  the  lnvlscld  limit,  we  use  the  Method  of 
Characteristics  (MOC)  procedure  to  determine  tne 
number  and  type  of  boundary  conditions.  This  theory 
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makes  It  clear  that  the  variable  at  the  boundary 
must  be  determined  from  a  combination  of  formal 
boundary  conditions  augmented  by  a  subset  of  the 
equations  of  motion.  This  rigorous  understanding  is 
lost  when  viscous  effects  are  present.  In  cases 
where  the  inertial  terms  dominate,  we  drop  normal 
diffusion  gradients  and  apply  the  HOC  procedures  as 
though  the  flow  were  Invlscld.  For  the  present 
problem,  these  procedures  corresponded  to  specifying 
the  stagnation  temperature,  the  stagnation  pressure, 
and  the  flow  angle  at  the  Inflow  boundary  and  tne 
back  pressure  at  the  outflow  boundary. 

For  low  Reynolds  numbers,  the  procedure  at  the 
downstream  boundary  was  changed  to  reflect  the 
dominance  of  the  viscous  terms,  but  the  lack  of  a 
rigorous  procedure  for  determining  the  type  of 
boundary  conditions  leads  to  some  ambiguity  as  to  how 
they  are  applied.  The  common  assumption  Is  that 
there  will  oe  one  additional  boundary  condition  for 
each  addition  order  of  derivative.  Thus,  the  four 
boundary  conditions  of  invlscld  flow  Increase  to 
seven  boundary  conditions  In  viscous  flow.  In  the 
present  calculations,  the  Invlscld  upstr.-am  boundary 
conditions  given  above  were  retained  (with  the 
argument  that  we  were  treating  slug  flow  at  the 
Inlet),  and  the  downstream  specification  of  pressure 
was  augmented  by  setting  these  second  derivatives  to 
zero, 
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so  that  four  boundary  conditions  were  specified  at 
the  outflow  boundary. 

RESULTS  AND  DISCUSSION 


*11  results  presented  are  based  on  the  direct 
solution  procedure.  The  comparison  of  viscous  and 
Invlscld  boundary  conditions  is  discussed  first. 
Intuitive  reasoning  would  suggest  that  the  InviaclC 
boundary  conditions  would  hold  at  very  high  Reynolds 
numbers  while  the  viscous  boundary  conditions  would 
provide  accurate  solutions  at  very  low  Reynolds 
numbers.  It  would  be  hoped  that  both  sets  of 
boundary  conditions  would  provide  similar  results  for 
a  range  of  Intermediate  Reynolds  numbers.  To 
ascertain  the  behavior  or  these  boundary  conditions, 
numerical  experiments  were  run  at  Reynolds  numbers  of 
5,  50  and  500.  In  the  absence  of  buoyancy,  both  sets 
of  boundary  conditions  gave  similar  convergence  and 
similar  solutions  (the  solutions  were 
Indistinguishable  on  a  plot)  at  all  three  Reynolds 
numbers.  In  the  presence  of  buoyancy,  all  attempts 
with  the  invlscld  boundary  conditions  diverged  at 
Reynolds  numbers  below  500  for  high  heat  addition, 
but  the  viscous  boundary  conditions  provided 
effective  convergence.  * s  Is  shown  later,  the  effect 
of  buoyancy  Is  almost  negligible  at  Reynolds  numbers 
above  500  and  becomes  Increasingly  dominant  at  lower 
Reynolds  numbers.  (Recall,  Reynolds  number  Is  being 
controlled  by  changing  tne  velocity.) 

Flowfleld  results  for  the  duct  flow  with  heat 
addition  area  given  In  Figs.  3”&.  Results  are 
presented  for  calculations  with  and  without  buoyancy, 
the  Prandtl  number  being  held  constant  at  0.7  for  all 
cases. 


Figures  3  and  *  present  the  temperature  and 
velocity  profiles  at  the  exit  for  the  non-buoyant 
case  at  Reynolds  number  of  5,  50  and  500.  In  tne 
absence  of  buoyancy  the  heat  source  has  a  passive 
behavior  that  only  causes  an  expansion  of  tne  gas. 
Low  velocities  abet  the  rapid  diffusion  of  heat  from 
the  centerline  to  the  wall,  so  that  the  temperature 


becomes  essentially  uniform  at  a  Reynolds  number  of  5. 
At  the  higher  Reynolds  numbers,  tne  reduced  diffusion 
rates  lead  to  steeper  temperature  gradients,  with  the 
heating  effect  being  local  lied  to  the  center  of  the 
flow.  The  velocity  profiles  on  Fig.  A  corroborate 
these  Inferences.  The  rapid  diffusion  at  the  low 
Reynolds  number  case  leads  to  a  parabolic,  fully 
developed  profile.  The  temperature  gradients  at 
higher  Reynolds  numbers  give  rise  to  a  fuller 
profile,  Indicating  a  flow  which  Is  still  developing. 

The  corresponding  profiles  with  buoyancy 
Included  are  shown  In  Figs.  5  and  6.  The  buoyancy 
effect  Is  negligible  at  Reynolds  number  500,  with  the 
temperature  and  velocity  profiles  being  almost 
Identical  to  the  non-buoyant  ease,  as  can  be  seen  by 
comparing  Figs.  3  and  5  with  Figs.  2  and  A.  The 
dramatic  effect  of  buoyancy  Is,  however,  evident  at 
the  two  lower  Reynolds  numbers.  The  center  of  the 
flew  is  accelerated  due  to  heating,  and  the 
velocities  at  the  centerline  In  dimensional  form 
become  comparable  In  magnitude  to  the  velocity  at 
Reynolds  number  500. 

The  sharp  decrease  In  the  rate  of  diffusion  at 
these  lower  Reynolds  numbers  Is  reflected  In  the 
modest  temperature  rise  at  the  wall.  In  contrast  to 
the  non-buoyant  case,  the  temperature  profile  with 
buoyancy  at  a  Reynolds  number  of  5  Is  now  highly 
non-uniform.  The  velocity  profiles  also  differ  quite 
strongly  from  the  fully  developed  profiles  In  the 
non-buoyant  flow.  A  fairly  wide  region  near  the  wall 
has  very  low  velocities  with  a  tendency  for  flow 
reversal. 


CONCLUSIONS 


Time  dependent  numerical  algorithms  were  tested 
for  low  Mach  number,  low  Reynolds  number  flows  with 
and  without  buoyancy.  Results  show  that  approximate 
factorization  techniques  do  not  converge  for  cases  of 
interest.  Stability  analyses  or  these  schemes 
indicate  that  the  reason  ror  the  impaired  performance 
is  the  hybrid  At2  error  term  that  arises  Trom  the 
product  of  viscous  and  lnvlsclo  derivatives. 

Stability  analyses  for  the  unfactored  algorithm 
Indicate  unconditional  stability  for  all  time  step 
slies,  and  suggest  that  very  large  At's  will  give 
maximum  convergence  when  direct  Inversion  procedures 
are  used.  The  robustness  of  the  direct  Inversion 
routine  reinforces  the  Inference  that  the  instability 
encountered  at  low  Reynolds  numbers  Is  because  or  tne 
approximate  factorization. 

Calculations  with  both  Invlscld  and  viscous 
boundary  conditions  at  the  outflow  boundary  show  a 
broad  region  of  overlap  In  the  ranges  of  application 
for  the  two  procedures  where  buoyancy  Is  absent.  In 
the  presence  of  buoyancy,  the  Invlscld  boundary 
conditions  failed  to  provide  solutions  In  Reynolds 
number  ranges  where  buoyant  effects  were  significant. 
The  reason  for  this  was  traced  to  a  tendency  for 
reentry  flow  near  the  wall  at  the  exit  plane. 


Flowfleld  solutions  show  that  buoyancy  begins  to 
dominate  at  a  Reynolds  number  of  about  500.  The  body 
force  due  to  buoyancy  causes  the  velocity  on  the 
centerline  to  Increase  as  compared  with  the 
non-buoyant  calculations  and  this  Increased  speed 
decreases  the  heat  flux  to  the  walls. 
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Fig.  3  Nondlmenslonal  Temperature  Profiles  at  Exit 
Plane  Without  Buoyancy. 

9  -  (T-Tq^/CThi-Tcl) 


Radial  Oistance,  Y 


Fig.  a  Nondlmenslonal  Velocity  Profiles  at  Exit 
Plane  Without  Buoyancy. 
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Fig.  5  Nondlmenslonal  Temperature  Profiles  at  Exit 
Plane  With  Buoyancy. 
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Fig.  6  Velocity  Profile  at  Exit  Plane  With  Buoyancy 
Effects  Incluoed. 
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